plot_ly() and ggplotly() functionsplot_geo()We will work with COVID data downloaded from the New York Times. The dataset consists of COVID-19 cases and deaths in each US state during the course of the COVID epidemic.
The objective of this lab is to explore relationships between cases, deaths, and population sizes of US states, and plot data to demonstrate this
cv_states_readin <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv") )
state_pops <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv"))
state_pops$abb <- state_pops$state
state_pops$state <- state_pops$state_name
state_pops$state_name <- NULL
cv_states <- merge(cv_states_readin, state_pops, by="state")
head, and tail of
the datadim(cv_states)
## [1] 39114 9
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 1 Alabama 2021-08-21 1 659750 12000 1 4887871 96.50939 AL
## 2 Alabama 2021-09-17 1 764839 13048 1 4887871 96.50939 AL
## 3 Alabama 2020-07-24 1 76005 1438 1 4887871 96.50939 AL
## 4 Alabama 2021-09-03 1 714860 12394 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
## 6 Alabama 2022-02-09 1 1254032 17452 1 4887871 96.50939 AL
tail(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 39109 Wyoming 2021-04-09 56 56873 701 56 577737 5.950611 WY
## 39110 Wyoming 2021-09-18 56 83958 918 56 577737 5.950611 WY
## 39111 Wyoming 2020-06-01 56 910 17 56 577737 5.950611 WY
## 39112 Wyoming 2021-04-20 56 57456 705 56 577737 5.950611 WY
## 39113 Wyoming 2020-04-11 56 343 0 56 577737 5.950611 WY
## 39114 Wyoming 2021-01-06 56 45890 464 56 577737 5.950611 WY
str(cv_states)
## 'data.frame': 39114 obs. of 9 variables:
## $ state : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ date : IDate, format: "2021-08-21" "2021-09-17" ...
## $ fips : int 1 1 1 1 1 1 1 1 1 1 ...
## $ cases : int 659750 764839 76005 714860 51 1254032 159169 932250 3953 104786 ...
## $ deaths : int 12000 13048 1438 12394 0 17452 2558 16503 114 1882 ...
## $ geo_id : int 1 1 1 1 1 1 1 1 1 1 ...
## $ population : int 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
## $ pop_density: num 96.5 96.5 96.5 96.5 96.5 ...
## $ abb : chr "AL" "AL" "AL" "AL" ...
cv_states$date <- as.Date(cv_states$date, format="%Y-%m-%d")
state_list <- unique(cv_states$state)
cv_states$state <- factor(cv_states$state, levels = state_list)
abb_list <- unique(cv_states$abb)
cv_states$abb <- factor(cv_states$abb, levels = abb_list)
cv_states = cv_states[order(cv_states$state, cv_states$date),]
str(cv_states)
## 'data.frame': 39114 obs. of 9 variables:
## $ state : Factor w/ 52 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ date : Date, format: "2020-03-13" "2020-03-14" ...
## $ fips : int 1 1 1 1 1 1 1 1 1 1 ...
## $ cases : int 6 12 23 29 39 51 78 106 131 157 ...
## $ deaths : int 0 0 0 0 0 0 0 0 0 0 ...
## $ geo_id : int 1 1 1 1 1 1 1 1 1 1 ...
## $ population : int 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
## $ pop_density: num 96.5 96.5 96.5 96.5 96.5 ...
## $ abb : Factor w/ 52 levels "AL","AK","AZ",..: 1 1 1 1 1 1 1 1 1 1 ...
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
## 380 Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
## 431 Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
## 33 Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
## 168 Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
tail(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 38547 Wyoming 2022-03-18 56 155907 1769 56 577737 5.950611 WY
## 39053 Wyoming 2022-03-19 56 155907 1769 56 577737 5.950611 WY
## 38906 Wyoming 2022-03-20 56 155907 1769 56 577737 5.950611 WY
## 38959 Wyoming 2022-03-21 56 155907 1769 56 577737 5.950611 WY
## 39020 Wyoming 2022-03-22 56 155988 1783 56 577737 5.950611 WY
## 38587 Wyoming 2022-03-23 56 155988 1783 56 577737 5.950611 WY
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
## 380 Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
## 431 Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
## 33 Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
## 168 Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
summary(cv_states)
## state date fips cases
## Washington : 793 Min. :2020-01-21 Min. : 1.00 Min. : 1
## Illinois : 790 1st Qu.:2020-09-05 1st Qu.:16.00 1st Qu.: 49810
## California : 789 Median :2021-03-12 Median :29.00 Median : 209010
## Arizona : 788 Mean :2021-03-12 Mean :29.78 Mean : 536429
## Massachusetts: 782 3rd Qu.:2021-09-16 3rd Qu.:44.00 3rd Qu.: 665336
## Wisconsin : 778 Max. :2022-03-23 Max. :72.00 Max. :9079164
## (Other) :34394
## deaths geo_id population pop_density
## Min. : 0.0 Min. : 1.00 Min. : 577737 Min. : 1.292
## 1st Qu.: 884.2 1st Qu.:16.00 1st Qu.: 1805832 1st Qu.: 43.659
## Median : 3525.5 Median :29.00 Median : 4468402 Median : 107.860
## Mean : 8812.6 Mean :29.78 Mean : 6419723 Mean : 422.717
## 3rd Qu.:10688.8 3rd Qu.:44.00 3rd Qu.: 7535591 3rd Qu.: 229.511
## Max. :88461.0 Max. :72.00 Max. :39557045 Max. :11490.120
## NA's :741
## abb
## WA : 793
## IL : 790
## CA : 789
## AZ : 788
## MA : 782
## WI : 778
## (Other):34394
min(cv_states$date)
## [1] "2020-01-21"
max(cv_states$date)
## [1] "2022-03-23"
new_cases and new_deaths and
correct outliersnew_cases, and new deaths,
new_deaths:
new_cases equal to the difference
between cases on date i and date i-1, starting on date i=2for (i in 1:length(state_list)) {
cv_subset = subset(cv_states, state == state_list[i])
cv_subset = cv_subset[order(cv_subset$date),]
# add starting level for new cases and deaths
cv_subset$new_cases = cv_subset$cases[1]
cv_subset$new_deaths = cv_subset$deaths[1]
for (j in 2:nrow(cv_subset)) {
cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j-1]
cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j-1]
}
cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases
cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths
}
# Focus on recent dates
cv_states <- cv_states %>% dplyr::filter(date >= "2021-07-01")
ggplotly for EDA: See if there are outliers or
values that don’t make sense for new_cases and
new_deaths. Which states and which dates have strange
values?p1 <- ggplot(cv_states, aes(x = date, y = new_cases, color = state)) +
geom_line()
ggplotly(p1)
There were -4678 recorded new cases on January 29, 2022 in Colorado, and -4397 recorded new cases in Pennsylvania on February 8, 2022. Neither of these values seem realistic, as we cannot have negative amounts of new cases.
p2 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +
geom_line()
ggplotly(p2)
There are a few instances in which a negative amount of people were recorded as dying on a particular day, such as -3770 deaths in Massachusetts on March 14, 2022, or -357 deaths in California on August 11, 2021.
Additionally, there are other outliers in the data where some states reported incredibly high numbers on a particular day compared to other days, such as the single-day spikes for Missouri and Tennessee.
Correct outliers: Set negative values for new_cases
or new_deaths to 0
Inspect data again interactively
cv_states$new_cases[cv_states$new_cases < 0] = 0
cv_states$new_deaths[cv_states$new_deaths < 0] = 0
p3 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +
geom_line() + geom_point(size = 0.5, alpha = 0.5)
ggplotly(p3)
Add population-normalized (by 100,000) variables for each
variable type (rounded to 1 decimal place). Make sure the variables you
calculate are in the correct format (numeric). You can use
the following variable names:
per100k = cases per 100,000 populationnewper100k= new cases per 100,000deathsper100k = deaths per 100,000newdeathsper100k = new deaths per 100,000Add a “naive CFR” variable representing
deaths / cases on each date for each state
Create a dataframe representing values on the most recent date,
cv_states_today
cv_states$per100k = as.numeric(format(round(cv_states$cases/(cv_states$population/100000),1),nsmall=1))
cv_states$newper100k = as.numeric(format(round(cv_states$new_cases/(cv_states$population/100000),1),nsmall=1))
cv_states$deathsper100k = as.numeric(format(round(cv_states$deaths/(cv_states$population/100000),1),nsmall=1))
cv_states$newdeathsper100k = as.numeric(format(round(cv_states$new_deaths/(cv_states$population/100000),1),nsmall=1))
cv_states = cv_states %>% mutate(naive_CFR = round((deaths*100/cases),2))
max_date <- max(cv_states$date)
cv_states_today = cv_states %>% filter(date==as.Date(max_date))
plot_ly()plot_ly() representing
pop_density vs. various variables (e.g. cases,
per100k, deaths, deathsper100k)
for each state on most recent date (cv_states_today)
cv_states_today %>%
plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
color = ~state, size = ~population, sizes = c(5, 70),
marker = list(sizemode = 'diameter', opacity = 0.5))
Washington D.C. is a massive outlier compared to the 50 states, so we will remove it and plot the remaining areas.
cv_states_today %>%
filter(state != "District of Columbia") %>%
plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
color = ~state, size = ~population, sizes = c(5, 70),
marker = list(sizemode = 'diameter', opacity = 0.5))
hovermode = "compare"cv_states_today %>%
filter(state != "District of Columbia") %>%
plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
color = ~state, size = ~population, sizes = c(5, 70),
marker = list(sizemode = 'diameter', opacity = 0.5),
hover_info = "text",
text = ~paste0("State: ", state,
"<br>Cases per 100k: ", per100k,
"<br>Deaths per 100k: ", deathsper100k)) %>%
layout(title = "Population-normalized cases per 100k",
yaxis = list(title = "Cases per 100k"),
xaxis = list(title = "Population Density"),
hovermode = "compare")
ggplotly()pop_density vs. newdeathsper100k
create a chart with the same variables using
ggplotly()pop_density
correlates with newdeathsper100k?p4 <- cv_states_today %>%
filter(state != "District of Columbia") %>%
ggplot(aes(x = pop_density, y = deathsper100k, colour = state, size = population)) +
geom_point()
ggplotly(p4)
Population density and new deaths per 100k seem to have some amount of correlation with one another, as the deaths per 100k people in a state seems to increase as the population density of the state increases. However, the relationship seems to be non-linear, so there are likely multiple other factors at play in the relationship between these variables.
naive_CFR for all states
over time using plot_ly()
naive_CFR for
the states that had an increase in September. How have they changed over
time?new_cases and new_deaths together in one plot.
Hint: use add_layer()
cv_states %>%
plot_ly(x = ~date, y = ~naive_CFR, color = ~state, type = "scatter", mode = "lines")
cv_states %>%
filter(state == "Florida") %>%
plot_ly(x = ~date, y = ~new_cases, type = "scatter", mode = "lines") %>%
add_lines(x = ~date, y = ~new_deaths)
In many of the states, deaths seemed to peak around November 2021, whereas the new cases seemed to peak around January 2022. This is likely just because of the different strains of COVID which were rampant at these times, as the strain which was more common in January 2022 (Omicron) was significantly less deadly than earlier strains of the virus.
Create a heatmap to visualize new_cases for each state
on each date greater than July 1st, 2021 - Start by mapping selected
features in the dataframe into a matrix using the tidyr
package function pivot_wider(), naming the rows and
columns, as done in the lecture notes - Use plot_ly() to
create a heatmap out of this matrix. Which states stand out?
cv_states_mat <- cv_states %>%
select(state, date, new_cases) %>%
filter(date > "2021-07-01")
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state,
values_from = new_cases))
cv_states_mat2 <- cv_states_mat2 %>%
column_to_rownames('date') %>%
as.matrix()
plot_ly(x = colnames(cv_states_mat2), y = rownames(cv_states_mat2),
z = ~cv_states_mat2, type = 'heatmap')
The heatmap above looks almost entirely purple, and the only extraordinarily high values for cases per day seem to occur in the states with higher populations (California, New York, Texas). This information isn’t that relevant on its own, as we’d expect states with high populations to have higher case counts.
new_cases for each state over time becomes more clear by
filtering to only look at dates every two weeks (check out the omicron
wave).filter_dates <- seq(as.Date("2021-11-01"), as.Date("2022-03-22"), by="2 weeks")
cv_states_mat <- cv_states %>% select(state, date, new_cases) %>% filter( date %in% filter_dates )
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = new_cases))
rownames(cv_states_mat2) <- cv_states_mat2$date
cv_states_mat2$date <- NULL
cv_states_mat2 <- as.matrix(cv_states_mat2)
#create heatmap
naive_CFR by state on
Devember 15, 2021naive_CFR by state
on most recent datepick.date = "2021-12-15"
# Extract the data for each state by its abbreviation
cv_per100 <- cv_states %>% filter(date==pick.date) %>% select(state, abb, newper100k, cases, deaths) # select data
cv_per100$state_name <- cv_per100$state
cv_per100$state <- cv_per100$abb
cv_per100$abb <- NULL
# Create hover text
cv_per100$hover <- with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))
# Set up mapping details
set_map_details <- list(
scope = 'usa',
projection = list(type = 'albers usa'),
showlakes = TRUE,
lakecolor = toRGB('steelblue')
)
# Make sure both maps are on the same color scale
shadeLimit <- 125
# Create the map
# Map for today's date
# Extract the data for each state by its abbreviation
cv_per100 <- cv_states_today %>% select(state, abb, newper100k, cases, deaths) # select data
cv_per100$state_name <- cv_per100$state
cv_per100$state <- cv_per100$abb
cv_per100$abb <- NULL
# Create hover text
cv_per100$hover <- with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))
# Set up mapping details
set_map_details <- list(
scope = 'usa',
projection = list(type = 'albers usa'),
showlakes = TRUE,
lakecolor = toRGB('steelblue')
)
# Create the map
Lab 10b questions 1-2, lab 11 questions 0-9 (only first half of q9). Upload html or pdf for both lab Rmd’s to quercus.